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ABSTRACT 



Aims. The extragalactic background light (EBL) in the ultraviolet to far-infrared wavelength region carries important information 
about galaxy and star formation history. Direct measurements are difficult, especially in the mid infrared region. We derive limits on 
the EBL density from the energy spectra of distant sources of very high energetic y-rays (VHE y-rays). 

Methods. The VHE y-rays are attenuated by the photons of the EBL via pair production, which leaves an imprint on the measured 
spectra from distant sources. So far, there are 14 detected extragalactic sources of VHE y-rays, 13 of which are TeV blazars. With 
physical assumptions about the intrinsic spectra of these sources, limits on the EBL can be derived. In this paper we present a new 
method of deriving constraints on the EBL. Here, we use only very basic assumptions about TeV blazar physics and no pre-defined 
EBL model, but instead a large number of generic shapes constructed from a grid in EBL density vs. wavelength. In our study we 
utilize spectral data from all known TeV blazars, making this the most complete study so far. 

Results. We derive limits on the EBL for three individual TeV blazar spectra (Mkn501, H 1426+428, 1ES 1101-232) and for all 
spectra combined. Combining the results from individual spectra leads to significantly stronger constraints over a wide wavelength 
range from the optical (~ 1 fim) to the far-infrared (~ 80 //m). The limits are only a factor of 2 to 3 above the absolute lower limits 
derived from source counts. In the mid-infrared our limits are the strongest constraints derived from TeV blazar spectra so far over an 
extended wavelength range. A high density of the EBL around 1 /mi, reported by direct detection experiments, can be excluded. 
Conclusions. Our results can be interpreted in two ways, (i) The EBL is almost resolved by source counts, leaving only a little room 
for additional components, such as the first stars, or (ii) the assumptions about the underlying physics are not valid, which would 
require substantial changes in the standard emission models of TeV blazars. 

Key words. BL Lacertae objects: general - diffuse radiation - galaxies: active - gamma rays: observations - infrared: general 
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The diffuse extragalactic background light (EBL) in the UV to 
far-IR wavelength regime carries important information about 
the galaxy and star formation history of the universe. The present 
EBL consists of the integrated electromagnetic radiation from 
all epochs, which is redshifted, corresponding to its formation 
epoch. In the energy density distribution, a two-peak structure 
is commonly expected: a first peak around 1 /mi produced by 
starlight and a second peak at ~100//m resulting from starlight 
that has been absorbed and reemitted by dust in galaxies. Other 
contributions, like emission from AGN and quasars, are ex- 
pected to produce no more t han 5 to 20% of the total EBL den- 
sity in the mid IR (see e.g. Mat ute et al.l 12006 and references 
therein). 

Solid lower limi ts for the EBL level have been derived from 
source counts (e.g . iMadau & Pozzettl 2000; Fa zio et all l2004t 
iFraver et al.l 12006; Dol e et al.l 120061) . Direct measurements of 
the EBL have proven to be a difficult task due to dominant 
fore grounds mainly fro m inter-planetary dust (zodiacal light) 
(e.g. lHauser et alJ ll998). and it is not expected that the sensi- 
tivity of measurements will greatly improve over the next few 
years. Several upper limits were reported from direct observa- 
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tions (e.g.lHauser et alj! 998) and from fluctuation analyses (e.g. 
iKashlinskv & Odenwaldll2000l) . In total, the collective limits on 
the EBL between the UV and far-IR confirm the expected two- 
peak structure, although the absolute level of the EBL density 
remains uncertain by a factor of 2 to 10 (see Fig.[TJ. In addition, 
several direct detections have also been reported, which do not 
contr adict the limits but lie sig nificantly above th e lower limits 
(see lHauser & Dwekll200ll and IKashlinskv! 120051 for recent re- 
views). In particular, a diffuse residual excess in the near IR (1 
to_4_um) was reported by the IRTS satellite dMatsumoto et al.l 
2005), which is significantly higher than the EBL density ex- 
pected from galaxy number counts. 

The reported excess lead to a controversial discussion about 
its origin. If the excess were extragalactic origin (i.e. associ- 
ated with the EBL), it might be attributed to emissions by the 
first stars in the history of the universe (Population III) and 
would make the EBL and its structure a unique probe of the 
epoch of Population III form ation and evolution ( Mapell Tet al.l 
12004 IKashlinskv et al.ll2005l) . Such a luminous Population III 
star generation, however, over- p redicts the number of Ly-g e mit- 
ters in ultra-deep field searches (Salva terra & Ferrarall2006l) and 
violates common assumptions on ba ryon consumpti on and star 
formati on r ates dDwek et al.| l2005al) . In addition, iDwek et all 
(20053) and lMatsumoto et all d2005l) point out that the NIR ex- 
cess could be attributed to zodiacal light. 
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Fig. 1. EBL measurements and limits. Upp e r lim- 

(120001) 
) (open 



its in the UV to optical: EdelsteinetaL 
(grey filled triangle - ). iMartin et alT l 1 99 1 



pink circle). iBrown e t al. (2000) (filled pink trian- 
gle). iMattilal d 19901) ( open g reen triangle), iToller! 
dl983l) / ILeinert et all d 19981) (open g reen square), 
iDube et all (1 19791) / ILeinert et all (1 1998b (open green 
diamond); Ten ta tive d e tectio n in the UV/optical: 
iBernstein et ail d2002l 120051) ( filled red circle); 
Lower limits from source counts: Madau & Pozzetti 
d2000b (open grey tri angles), iFazioetall d2004l) 
(open blue triangles ) , lElbaz et al.1 d2002l) (green 
cross) iMetcalfe et all d2003l) (red x).lPapovich et al.1 
d2004l) (filled re d triangle), iDole et alj d2006b (filled 
pink triangles), iFraver et al.l d2006j) (open red tri- 
angle) ; Detections in the nea r IR: iDwek & Arendtl 
(1998) (open pink cross), Gor iian et all d200C' 
(filled brown circle), Wrig ht & Reesd 
(open blue squares), Cambresv et al.l 
(filled brown squares), Matsumoto et al.l 
(small open grey circles), iLevenson et al.l 
(open blue ci rcles); Uppe r limit s from direct 
measurements: lHauser et al.l (11998b (filled green 
triangles), IDwek & Arendtl d 1998b (filled pink 
triangles), Lag ache & Pugetl d2000l) (filled blue 
tria ngles); Upper limits fro m fluctuation analy- 

blue 




sis: Kashlinsky et al. (1996) (filled blue circles), 
iKashlinskv & Odenwaldl d2000b (filled pink circles); 
Lower hm i ts fro m stacking analysis in the far-IR: 
IDole et all ( 2006) (blue triangles); Detections in the 
far-IR: | Hauseret aD (19 98) (filled green squares), 
lLagache & Pugetl d2000b (ten tative, filled blue 
square). iFinkbeiner et al.l (2000) (tentative, open red 
diamonds). 



A different app roach of deriving co nstraints on the EBL (la- 
beled "clever" bv lDwek & Slavir3ll994b became available with 
the detec tion of very high -energy (VHE) y-rays from distant 
sources dPunch et al.lll992b . These VHE y-rays are attenuated 
via pair production with low-energy photons from the EBL 
dGould & Schrederlll967h . It was soon realized that the mea- 
sured spectra can be used to test the transparency of the uni- 
verse to VHE y-rays and thus t o derive constraints on the EBL 
density dFazio & Steckeriri970b . With reasonable assumptions 
about the intrinsic spectrum emitted at the source, limits on the 
EBL density can be inferred. In a pioneering work conducted by 
IStecker et all (1996), first limits on the EBL were derived. The 
method is, however, not straightforward. It is, in principle, not 
possible to distinguish between source-inherent effects (such as 
absorption in the source, highest particle energies, magnetic field 
strength, etc.) and an imprint of the EBL on a measured VHE 
spectrum. The emission processes in the detected extragalactic 
VHE y-ray emitter^] are far from being fully understood, which 
makes it more difficult to use robust assumptions on the shape 
of the intrinsic spectrum. Different EBL shapes can lead to the 
same attenuation imprint in the measured spectra, thus a recon- 
structed attenuation imprint cannot uniquely be identified with 
one specific EBL shape. A further uncertainty arises due to the 
unknown evolution of the EBL with time, which becomes more 
important for distant sources. 



1 So far, all but one detected extragalactic VHE y-ray emitters belong 
to the class of TeV blazars. 



Nevertheless, measured VHE y-ray spectra can provide ro- 
bust upper limits on the density and spectral distribution of the 
EBL, if conservative assumptions about the emission mecha- 
nisms of y-rays are considered. A review of the various efforts 
to detect the EBL or to derive upper limits via the observa- 
tion of extragala ctic VHE y-ray s ource s up to the year 2001 
can be found in Ha user & Dwekl d200lb . The measured spec- 
trum of the TeV blazar H 1426+428 at a redshift of z = 0.129 
dAharonian et alj|2002bb led to a first t entative detection of an 
imprin t of the EBL in a VHE spectrum ( Aharon ian et alj |2002b, 
2003c). Using the H 1426+428 spectrum, together with the spec- 
tra of pr eviously detected TeV bl a zars, limits on the E BL were 
derived dCostamante et alj 120031; iKneiske et ail 12004b . Later, 
IDwek & Krennrichl d2005l) utilized a large sample of TeV blazar 
spectra and solid statistical methods to test a set of EBL shapes 
on their physical feasibility. The EBL shapes were considered 
forbidden, when, under the most conservative assumptions, the 
intrinsic spectra showed a significant exponential rise at high 
energies. Using the VHE spectra of the T eV blazars Mkn421 , 
Mkn501, H 1426+428, and PKS 2155-304. lDwek et alJd2005bb 
argued that the claimed NIR excess is very unlikely to be of 
extragalactic origin and that an EBL density on the level of 
the source counts gives a good representation of the intrinsic 
spectra. Recently, the detection of the two distant TeV blazars 
H 2356-309 and 1ES 1101-232 by the H.E.S.S. experiment has 
been used to derive strong li mits on the EBL density around 
2 /mi dAharonian et al. 2006a). The method is based on the hy- 
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pothesis that the intrinsic TeV blazar spectrum cannot be harder 
than a theoretical lim it and that the spect rum of the EBL fol- 
lows a certain shape dPrimack et al.|[l999l) . To derive limits on 
the EBL density, the shape is scaled until the de-absorbed TeV 
blazar spectrum meets the exclusion criteria^ The derived upper 
limits are only a factor of two above the lower limits from the 
integrated light of resolved galaxies. 

In the past 2-3 years, many new TeV blazars have been 
discovered and the established ones were remeasured with 
higher sensitivity with the new Imaging Atmospheric Cherenkov 
Telescopes (IACT) such as H.E.S.S and MAGIC. It is therefore 
of high interest whether the new measurements agree with the 
previous limits on the EBL. Furthermore, to derive limits on 
the EBL from the data, it is important to scrutinize all avail- 
able data together to obtain a consistent picture and to maximize 
the constraints. In this paper we use, therefore, spectra from all 
TeV blazars to derive upper limits on the EBL density in a wide 
wavelength range from the optical to far-infrared. Moreover, a 
common criticism of the EBL limits previously derived with this 
method is that the limits are obtained by assuming a certain EBL 
model and e.g. scaling it, or by exploring just a few model pa- 
rameters. Since EBL models are complex and different models 
do not agree in details, the derived limits become very model- 
dependent. In order to avoid this dependency, we developed a 
novel technique of describing the EBL number density by spline 
functions, which allows us to test a large number of hypothetical 
EBL shapes. Our aims are to 

1 . Provide limits on the EBL density, which do not rely on a 
predefined shape or model, but rather allow for any shape 
compatible with the current limits from direct measurements 
and model predications. 

2. Treat all TeV blazar spectra in a consistent way, using sim- 
ple and generic assumptions about the intrinsic VHE y-ray 
spectra and a statistical approach to find exclusion criteria. 

3. Use spectral data from all detected TeV blazars to (a) derive 
upper limits on the EBL density from the individual spectra 
and then (b) combine these results into a single robust limit 
on the EBL density for a wide wavelength range. 

The paper is organized as follows: In Section|2]we describe 
our method of using splines and a grid in EBL density vs. wave- 
length to construct different shapes, in Section [3] we introduce 
the TeV blazar spectra, and in Section|4]we describe the criteria 
we impose on these spectra to derive limits on the density of the 
EBL. In Section[5]results for individual sources and in Section|6] 
the combined results are discussed. Last but not least, we give a 
conclusion and outlook in Section [7] 

Throughout the paper we adopt a Hubble constant of Hq = 
72kms~ 1 Mpc~ 1 and a flat universe cosmology with a matter 
density normalized to the critical density of Q. m = 0.3 and 
Q A = 0.7. 
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Fig. 2. Examples for spline shapes resulting from the grid lay- 
out, overlaid on the grid points (red filled circles). The dashed 
black line illustrates the thinnest structure that can be achieved 
with our grid setup. A Planck spectrum (blue line) is given for 
comparison. Further examples for shapes are given as solid black 
lines. 
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Fig. 3. Top: Grid points (red filled circles) and minimum and 
maximum shape of the scan. Bottom: Minimum and maximum 
shape overlaid on the EBL measurements from Fig.Q] 



2. Grid scan of the EBL with splines 

TeV y-rays traversing the extragalactic radiation fields are ab- 
sorbed via pair production with the low-energy photon s of the 
EBL: yxev Tebl — * <? + e~. Details on the cross section dHeitlerl 
1960) and the exact calc ulations can be found elsewhere (see 
Dwek & Krennrichll2005l for an overview). Here, we focus on 

2 A similar technique (scaling of a certain shape until the de- 
ab s rbed spec t rum re aches an exclusion criterion) was previously used 
by iGuv et al.l d200Qh to derive limits on the EBL from the Mkn 501 
spectrum measured by the CAT experiment during a flare in 1997. 



the new technique using splines (as introduced below) to be able 
to examine as many as several million possible EBL shapes. 

To calculate the optical depth t 7 for a TeV-y ray of energy 
E y emitted at a redshift of z for one realization/shape of the EBL, 
one needs to solve a three-fold integral over the distance I, the 
interaction angles between the two photons fi - cos 9 and the 
number density of the EBL photons n e \ 

Ty(E y , Z) = (1) 

JT d£(z) £ dju £ de' n e (e',z) cr yy (E' y , e', fx) 
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where <x rr is the pair production cross section, e the energy of 
the EBL photon, and E' y and e' refer to redshifted energy values. 
For the large number of shapes (~8 million) we will analyze, 
a full numerical integration would require an extensive amount 
of computing power. To avoid this problem, we parametrize the 
EBL number density at z — as a spline 



i=0 



with 



s;,o(<0 = 



1 if e, < e < e, + i and e, < e i+ \ otherwise 


e ~ £ i , \ e i+p+l ~ e , . 

Si,p-i(e) + s i+Up -i(e) . 



(2) 



(3) 



Here /? is the order of the spline and is set to p — 3 throughout 
the paper, k is the number of supporting points, e\, q are the 
positions of the supporting points of the curve, and w\, Wk are 
weights controlling the shape of the curve. 

A way to visualize this spline is to add a number of Gaussian- 
like base functions to lead to a smooth curve, whereby the 
base functions can be weighted (w,) to achieve different over- 
all shapes. These base functions are characterized by a position 
of the center and a certain width, which depends on the order p 
and the spacing of the supporting points. 

By inserting Eq. (O into Eq. (Q~|i and then swapping the inte- 
gration and the summation, one obtains an expression for the op- 
tical depth, where the integral can be solved easily for a certain 
redshift and set of supporting points of the spline. The optical 
depth can then be calculated by a simple summation, where the 
shape of the EBL is determined by the choice of a set of weights. 

To estimate the numerical uncertainties in the integration of 
the base functions, the result for r v from the summation is com- dN/dE^ = dN/dE b S x exp[r r (£' r , z)], 



scales, since they would be smoothed out by redshift. Of course 
the upper limits depend on the thickness of the structures mod- 
eled, and the choice of minimum thickness has to be physically 
motivated, as it is in our case. 

The weights w, of the spline (y-axis of the grid) are var- 
ied in 12 equidistant steps in log 1() (energy density) from 0.1 to 
100 nW irT 2 sr~' . Since we apply the grid point positions directly 
as weights w,, the resulting curve does not directly go through 
the grid points. So for all limits, etc., the actual curve has to 
be considered and not the grid point positions. Our use of the 
grid point positions as weights also results in several shapes ly- 
ing between two grid points. This is illustrated in Fig.[2]for two 
arbitrarily chosen grid points between 10 and 20 pm. 

For the scan, a subset of grid points is selected such that all 
resulting shapes are within the limits given by the galaxy counts 
on the lower end and by the limits from the direct measure- 
ments and the fluctuation analyses on the upper end. In the 20 
to lOOjurn regime, we choose a somewhat looser shape, given 
the wider spread of the (tentative) measurements (Fig. [3]). By it- 
erating over all grid points, we obtain 8 064 000 different EBL 
shapes, which will be examined. 

The generic shapes used in this analysis have no redshift de- 
pendency, thus evolution of the EBL cannot be taken into ac- 
count. We assume a constant photon number density only ex- 
panding and shifting in wavelength with expansion of the uni- 
verse. Neglecting the evolution of the EBL is a valid assumption 
for nearby sources (e.g. Mkn 501), but for distant sources this 
can result in an additional error in the order of 10 to 20% depend- 
ing on the evolution model and redshift (see e.g. A haronian et ail 
l2006al) . The possible systematic error arising from this assump- 
tion will be discussed further in Sect. [7] 

For a given optical depth t 7 (E 7 , z) one can determine the 
intrinsic differential energy spectrum using 



(4) 



pared to the results from a full integration over the actual EBL 
shapes for several shapes and spectra used in this paper. The de- 
viation for almost all settings is found to be less then 0.5%. Even 
for extreme cases (high redshift, high EBL density), the devia- 
tion is always less than 2%. These small deviations arise from 
inaccuracies in the numerical integration of the base functions. 

The spline parameterization is used to construct a set of EBL 
shapes using a grid in EBL energy density vs. wavelength. The 
x-positions of the grid points (wavelength) are used as positions 
for the supporting points e, of the spline. The y-positions of the 
grid points (energy density) are used as weights w,. 

For the supporting points (x-axis of the grid) e,, we use 16 
equidistant points in log 10 (/l) (A = hc/e) from 0.1 to 1000 pm. 
The number and distance of the supporting points, together with 
the order of the spline, determines the width of the structures, 
which can be described with the spline. Given the grid setup and 
the spline of order p = 3, the thinnest peak or dip that can be 
modeled is about three grid points wide (Fig. [2j. This minimum 
width is similar to a Planck spectrum from a black body radia- 
tion (Fig. [2j. It is expected that the EBL originates mainly from 
an overlap of the redshifted spectra of single stars, for which a 
Planck spectrum is a good first-order approximation, and dust 
reemission. Thus, an EBL structure thinner than a Planck spec- 
trum is unlikely. 

Noteworthy, any extremely sharp and strong cut-offs or 
bumps can only be described in an approximate way. Such 
features can arise from the Lyman-a drop-off of massive 
Popu lation III stars in certain models (e.g. ISalvaterra & Ferraral 
2003), but they are generally not expected on larger wavelength 



where dN/dE i, s is the observed spectrum. The intrinsic spec- 
trum is then tested on its physical feasibility as described in 
Section [4] This process is repeated for all VHE y-ray sources 
(described in the next section) and for all 8 064 000 EBL shapes. 



3. TeV blazar sample 

Since the detection of the fi rst extragalactic VHE y-ray source 
in 1992 dPunch et al.l fl992). a wealth of new data from this 
source type has become available. Up to now all detected VHE 
y-ray sources belong to the class of active galactic nuclei (AGN) 
and, with the only one exce ption of the radio galaxy M87 
(Aharonianet al. 2003b, 2006d), to the subclass of TeV blazars 



Horan & Weekes 2004). AGNs are known to be highly variable 
sources, with variations of the absolute flux levels in the TeV 
range by more than an order of magnitude and on time scales a s 
short as 15 min (e.g. lGaidos et alJl996tlAharonian et alj20 02a). 
Changes in the spectral shapes have also be en observed (e.g. 
iKrennrich et~aT]|2002t lAharonian etai1l2002al) . 

In this study we utilize spectral data obtained during the 
past seven years by four different experimental groups operating 
ground-based imaging atmospheric Cherenkov telescopes 
(IACT s): Whi pple dFinlev & The VERITAS Collaboration! 
120011). H EGRA dPaum et al J 19971) . H.E.S.S. dHintonll2004l) and 
MAGIC dCortina et alj|2005l) . We select at least one spectrum 
for every extragalactic source with known redshift. If there 
is more than one measurement with a comparable energy 
range, we take the spectrum with the better statistics and the 
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Source 


Redshift 


Experiment 


Energy range 


Slope 


Cut-off energy 


Reference 








(TeV) 


T ± cr s , ± cr sy 


(TeV) 




Mkn421 


0.030 


MAGIC 


0.10 


-3.0 


2.20 ± 0.08 ± 0.20 


1.44 ±0.28 


Albert et al. (2007b) 


Mkn421 


0.030 


HEGRA 


0.70- 


- 18.0 


2.19 ±0.02 + 0.20 


3.6 + 0.4 - 0.3 


Aharonian et al. (2002a) 


Mkn421 


0.030 


Whipple 


0.35- 


-0.90 


2.31 ±0.04 ±0.05 




Krennrich et al. (2002) 


Mkn501 


0.034 


HEGRA 


0.50- 


-22.0 


1.92 ±0.03 ±0.20 


6.2 ± 0.4 


Aharonian et al. (1999) 


1ES 2344+514 


0.044 


Whipple 


0.80- 


- 11.0 


2.54 ±0.17 ±0.07 





Schroedter et al. (2005) 


Mknl80 


0.045 


MAGIC 


0.14 


- 1.5 


3.25 ± 0.66 ± 0.20 




Albert et al. (2006b) 


1ES 1959+650 


0.047 


HEGRA 


1.5- 


13.0 


2.83 ±0.14 ±0.08 




Aharonian et al. (2003a) 


PKS 2005-489 


0.071 


H.E.S.S. 


0.20 


-2.5 


4.0 ±0.4 (±0.2) 




Aharonian et al. (2005a) 


PKS 2155-304 


0.116 


H.E.S.S. 


0.20 


-3.5 


3.37 ±0.07 ±0.10 




Aharonian et al. (2005b) 


H 1426+428 


0.129 


HEGRA 


0.70- 


- 12.0 


2.6 ±0.6 ±0.1 




Aharonian et al. (2003c) 


H 2356-309 


0.165 


H.E.S.S. 


0.16 


- 1.0 


3.06 ±0.21 ±0.10 




Aharonian et al. (2006b) 


1ES 1218+304 


0.182 


MAGIC 


0.08 


-0.7 


3.0 ± 0.4 ± 0.6 




Albert et al. (2006a) 


1ES 1101-232 


0.186 


H.E.S.S. 


0.16 


-3.3 


2.88 ±0.14 ±0.1 




Aharonian et al. (2006a) 
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1ES1101-232, H.E.S.S., x 1.7x1 (T 
1ES1 21 8+304, MAG IC,x 1.4x1 26 
H2356-309, H.E.S.S., x 2.0x1 25 
H1 426+428, HEGRA, x 2.8x1 21 
PKS 2155-304, H.E.S.S., x 7.8x1 2 



PKS 2005-489, H.E.S.S., x 1 .8x10' 



1 ES1 959+650, HEGRA, x 1 .7x1 1 ' 
Mkn 180, MAGIC, x 4.4x1 16 
1ES2344+514, Whipple, x 2.2x1 1 
Mkn 501, HEGRA, x 3.1x1 11 



-6 — Mkn 421, Whipple, x 3.6x10 



Mkn 421, HEGRA, x 1.1x1 1 
Mkn 421, MAGIC, x 1.5x1 1< 



Fig. 4. TeV blazar sample chosen for this study. From all detected TeV blazars, at least one spectrum is shown with the exception of 
PG 1553+1 13 (unknown redshift). The spectra are multiplied by E 2 to emphasize spectral differences and are spread out along the 
Y-axis and ordered in redshift to avoid cluttering of the plot (the corresponding scaling factors are given in the legend of the plot). 



harder spectrum (expected to give stronger constraints, see 
Sect. H]i. If different measurements of one source cover different 
energy ranges, we include both spectra as independent tests. 
Another possibility would be to combine the measurements 
from different e xperiments, as was don e before in the case 
of H 1426+428 dAharonian et alj 12003d) . However, since the 
sources are known to be variable in the flux level and spectral 
shape, a combination of non-simultaneous data is not trivial. 
In addition, in a conservative approach one has to consider the 
systematic errors reported by the individual experiments, which 
are around 20%. Our method is quite sensitive to the errors of 
the flux, and this additional error would weaken our results. We 
therefore do not use combined spectra. The selected TeV blazar 
spectra are summarized in Fig. [4] and Table [TJ In the case the 
measured spectrum can be described well by a simple power law 
dN/dE oc E~ r , only the slope F is quoted. Otherwise, the cut-off 
energy E co ff according to dN/dE oc E~ T exp(-£ , /£ , cof f) is also 



quoted. If no systematic error in the photon index is given in the 
corresponding paper, we use a value of 0.2 (values in brackets). 
In Fig. |4] the spectra are multiplied by E 2 to emphasize spectral 
differences and are spread out along the Y-axis and ordered in 
redshift to avoid cluttering the plot. 

We did not include data from the radio gala xy M 87, which 
was rec ently confirmed as a VHE y-ray emitter dAharonian et alj 
l2006dl) .M87 is a nearby source (z = 0.00436), so even for high 
EBL densities the attenuation is weak and would only be no- 
ticeable at high energies ~30 TeV. Nevertheless, M 87 has a hard 
spectrum with a photon index of F ~ 2.2 currently measured up 
to 20 TeV, so further observations that extend the energy range 
to even higher energies could make M 87 an interesting target 
for EBL studies as well. 

The cross section o" ry of the pair production process has a 
distinct peak close to the threshold. The maximum attenuation 
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of VHE photons of energy E y occurs by interaction with low- 
energy photons with a wavelength 

/T ax (jum) ~ 1.24 £ r (TeV). (5) 

Since the selected spectra cover an energy region from 100 GeV 
up to more than 20 TeV, the EBL wavelength range for the ab- 
sorption of VHE y-rays spans from UV (~0.1 yum) to the mid 
IR (~30ji/m). This region is of particular cosmological interest 
since it might contain a signature of Population III stars. 

4. Exclusion criteria for the EBL shapes 

We aim to construct an upper limit on the EBL density using 
the TeV blazar sample described in Sect. [3] In order to achieve 
this, we examine every EBL shape (as introduced in Sect. [2]) to 
see whether the intrinsic TeV blazar spectra, which result from 
correcting the measured spectra for the corresponding optical 
depths, are physically feasible. EBL shapes are considered to 
be allowed if the intrinsic spectra of all tested TeV blazars are 
feasible. As an upper limit we define the upper envelope of all 
allowed shapes. It is constraining in the wavelength range where 
it lies below the maximum shape of the scan. Otherwise no upper 
limit is quoted. 

There are different ways to examine a TeV blazar spectrum 
upon its feasibility. In this paper we follow very general argu- 
ments arising from the shock acceleration scenario of relativis- 
tic particles. In this well-accepted view, electrons are Fermi- 
accelerated with a resulting power-law spectrum of dN/dE ~ 
E~ a , with a slope a of about 2. Due to a faster cooling of high 
energy electrons compared to lower energies, the slope a can 
be steeper than 2, but there is no simple theoretical possibility 
to produce an electron spectrum with a harder spectrum. Given 
an electron spectrum, one can calculate the slope of the syn- 
chrotron energy spectrum: it is dN/dE ~ E~ r with a photon 
index F = 2 p- = 1.5. The energy spectrum of inverse-Compton 
photons, independent of the origin of the target photons, has ap- 
proximately the same photon index as the synchrotron energy 
spectrum if the scattering occurs in the Thomson regime. In the 
case of the Klein-Nishima regime, the index is even larger. Thus, 
the photon index of the energy spectrum of VHE photons orig- 
inating from an IC scattering is F = 1 .5 or larger. In the case 
of a hadronic origin, their spectrum is more complicated. Yet, 
the resulting VHE photons originate from pion decays leading 
to a photon spectrum with F = 2, which is softer than a pos- 
sible photon spectrum with a leptonic origin. In conclusion, we 
assume that the photon index F of the intrinsic TeV blazar spec- 
trum is 1.5 or larger. These ar gume nts were rec e ntly us ed in 
lAharonian eTail (I2003dl2006allcf) and lAlbert et all d2007al) . 

However, a possibility of obtaining even harder photon spec- 
tra is not fully excluded. For instance, thou gh contrary to the 
wide acceptance of these general arguments, Katarzvhski et al.l 
(2006) argue that synchrotron emission, as well as IC scattering 
of relativistic electrons, does not necessarily occur close to the 
region of electron acceleration. If so, the electron spectrum can 
become truncated due to propagation effects; i.e. the minimum 
energy of electrons can be as high as several GeV. In an extreme 
case, we deal with a monoenergetic spectrum of VHE relativistic 
electrons. Then, and this is the most extreme case, a resulting IC 
photon index can be as small as F = 2/3. We use this limit in 
the present paper, in addition to the standard limit of T = 1 .5, to 
demonstrate the strength of our method. 

In the simplest models, it is assumed that VHE photons origi- 
nate from a single compact region (so-called one-zone scenario). 



This result in a smooth, convex spectral energy-distribution with 
two peaks at certain energies in a log(vF(v))vs.log(v) represen- 
tation. Yet there are no obvious arguments against scenarios in- 
cluding some kind of multizone models, which can be naturally 
realized in the jets of TeV blazars. Then, the measured spectrum 
of VHE y-ray emitting sources will be a superposition of several 
one-zone emission regions. So far, there has been no indication 
of this in the measured spectra, however, the attenuation of VHE 
y-rays by EBL photons could hide such substructures. 

In addition to the constraints on the photon index of 
the intrinsic TeV blazar spectra, we also argue that a super- 
exponentially rising energy spectrum with an increasing energy 
is not reali stic. The so-called p i le-up a t high energies was first 
noticed by iProtheroe & Meverl ((2000) for the Mkn 501 spec- 
trum, and early att empts to avoid it invoked violation o f the 
Lorentz invariance (Kifune 1999; Proth eroe & Meyer! 2000). On 
the other hand, IStecker & Glashowl (1200 ll) argue that the same 
Mkn 501 spectrum data can be used to place severe constraints 
on the Lorentz-invariance breaking parameter. Another possi- 
bility for explaining pile-ups would require an ultrarelativistic 
j et with a very high bul k-motion Lorentz factor Fo > 3 x 10 7 
( Aharonia n et al.ll200 2c). The pile-ups, however, seem to arise 
at different energies for different sources, which is not expected 
in these models. Moreover, the pile-ups can be easily avoided by 
choosing a sufficiently low level of the EBL density. 

Based on the arguments above and including possible mul- 
tizone emission scenarios, we assume that at least one of the 
following smooth, analytical functions can describe the intrin- 
sic spectrum satisfactory well: a simple power law (PL), a bro- 
ken power law with a transition region (BPL), a broken power 
law with a transition region and a super-exponential pile-up 
(BPLSE), a double broken power law with two transition regions 
(DBPL), or a double broken power law with two transition re- 
gions and a super-exponential pile-up (DBPLSE). The functions 
are summarized in Table [2] 

In our approach of examining the intrinsic spectra of TeV 
blazars, we adopt a general assumption that one of the smooth 
functions above describes the intrinsic spectrum satisfactorily. 
In case we fail to find such a function, we abandon any further 
examination of the fit parameters, and a given EBL shape (real- 
ization) is not excluded. It is noteworthy that less than 0.06% of 
all intrinsic spectra of the TeV blazars could not be fitted well by 
the chosen functional forms, which are described below. 

In order to determine a good analytic description, we fit the 
intrinsic spectrum with the functions listed above, starting with 
the simplest one (PL). As "good description" we take a fit with 
a chance probability Ppit > 5% based on its^- 2 value. The deter- 
mined fit parameters are examined to be physically meaningful 
as described below. If the fit has a probability Ppit < 5%, we do 
not consider the function and examine the next function from the 
list (Table |2). Also if the fit has a chance probability Ppit > 5%, 
but the determined parameters do not lead to an exclusion of the 
corresponding EBL shape, we take the next function from the 
list. The reason is that a function with more free parameters can 
lead to a significantly better fit. In order to make sure that a more 
complicated function indeed describes the intrinsic s hape better 
than a simpler one, we use the likelihood ratio test (Eadi e et ail 
1988, Appendix A). We require at least a 95% probability to pre- 
fer a function with a higher number of free parameters. Only if 
the function is preferred, we examine its fitted parameters and 
their errors. Spectra with n data points are only fitted with func- 
tions with up to n - 1 parameters. In case of spectra with low 
statistics, we fix the softness of the break (/, see Table |2]i be- 
tween two spectral indices to an a priori chosen value of / = 4 
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Fig. 5. Limits on the EBL density from Mkn 501 (F max > 1.5). 
Grey solid curves are the minimum and maximum shapes of the 
scan; grey dashed curves are all the highest allowed EBL shapes 
for Mkn 501; the thick colored curve is the corresponding enve- 
lope shape. The grey shaded area marks the wavelength region, 
in which the envelope shape is constraining the EBL density. 

to allow for a fit by a higher-order function (e.g. in case of the 
BPL for H 1426+428). 

We finally define the following criteria to exclude an as- 
sumed EBL shape: 

- At least one of the determined photon indices from the best 
hypothesis is outside of the allowed range, i.e. F' max < r llmlt , 
where r limit = 1.5 or 2/3, and P- max = F + <x r , + o- sy with F 
and <x r the fitted photon index and its error, respectively, and 
o" sy as the systematic error on the spectral slope (as given for 
the corresponding measurement). 

- The best fit is obtained by one of the two shapes with an 
exponential pile-up (BPLSE or DBPLSE). 

A single confidence level for the derived upper limits on the 
EBL density cannot be quoted easily: For the test on the photon 
index we use a 1 cr confidence level (as defined for two-sided dis- 
tributions, including systematic errors). For the likelihood ratio 
test, we use a 2 cr level (as defined for one-sided distributions). 
Thus, the confidence level for the upper limit ranges from 68% 
to 95%. 

As discussed in this section, the theoretical expectations on 
the smallest possible photon index r llmlt have a given spread. 
Thus, we perform two EBL scans, assuming a limit of F lmlt = 
1.5 for the first realistic case and for the second extreme case 
a limit of F lmlt = 2/3. From here on, we refer to realistic and 
extreme scan accordingly. 

Since our exclusion criteria consider the effect of an EBL 
shape solely on the hardness of the spectra of extragalactic VHE 
y-ray sources, some shapes might be considered viable in our 
scan, even though other criteria for exclusion could be found 
(e.g. exclusion based on spectral shape of the EBL, related to 
the energy density of st ar and dust emission as discussed in 
iDwek & Krennrichll2005h . In this respect our approach is con- 
servative. 
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Fig. 6. Mkn 501 individual results (F max > 1.5). Upper Panel: 
Two examples for the highest allowed shapes (A - solid line; B - 
dashed line) overlayed on EBL measurements and the minimum 
and maximum shapes of the scan (grey). Middle Panel: Intrinsic 
spectrum for shape A and fit functions: BPL (solid black line), 
BPLSE (solid red line), DBPL (dashed black line), and DBPLSE 
(dashed red line). The fit parameters can be found in Table [3] 
The BPL and BPLSE are not good fits; the DBPL is a good fit 
and its parameters are within the allowed range; the DBPLSE is 
not prefered over the BPL and the DBPL. Lower Panel: Intrinsic 
spectrum for shape B and the same fit functions as in the Middle 
Panel. The fit parameters can be found in Table[3] All four func- 
tions are good fits. No function is preferred over the BPL. 



5. Results for individual spectra 

To present the results for individual spectra in a compact and 
non-repetitive way, we sort our preselected spectra into three 
categories. For each category we select one prototype spectrum 
(following similar criteria to those in Sect. 0), for which we give 
results. The categories considered are: 



"Nearby and well-measured". As prototype spectrum we select 
the Mkn 501 spectrum (z = 0.034) recorded by HEGRA dur- 
ing a major TeV flare in 1997. The relatively hard spectrum 
provides good statistics from 800 GeV up to 25 TeV. Other 
sources in this category are Mkn421 and 1ES 1959+650. 
1ES 2344+514 and Mkn 180 are at a comparable distance, 
but the spectra do not have such high statistics. PKS 2005- 
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Table 3. Fit results for the Mkn 501 spectrum using the functions 
from Fig. [6] 
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Fig. 7. Limits on the EBL density from H1426+428 (F max > 
1.5). Grey solid curves are the minimum and maximum shapes 
of the scan; grey dashed curves are all the highest allowed EBL 
shapes for H1426+428; the thick colored curve is the corre- 
sponding envelope shape. The grey shaded area marks the wave- 
length region, in which the envelope shape is constraining the 
EBL density. 



489 lies somewhere between this and the following category. 
Spectra in this category mainly provide limits in the FIR due 
to a pile-up at high energies. 
"Intermediate distance, wide energy range ". The prototype for 
this category is the spectrum from H 1426+428 at a red- 
shift of z = 0.129, with energies ranging from 700GeV up 
to 12 TeV. PKS 2155-304 is at a comparable distance, and 
its measured spectrum has a better statistic, but the highest 
energy point is only at 2.5 TeV; and the spectrum is much 
softer. 

"Distant source, hard spectrum ". The most distant TeV blazar 
discovered so far with a published energy spectrum is 
1ES 1101-232 at a redshift of z = 0.186. Its spectrum is 
hard and ranges from 160 GeV to 3.3 TeV. H 2356-309 and 
1ES 1218+304 are at similar distances, but the statistics 
are not as good and/or the spectrum is softer. This makes 
1ES 1 101-232 the natural choice for a prototype spectrum in 
this category. 

As the limit on the EBL for the individual spectra (and for 
the combined results in Sect. [6] as well), we define the enve- 
lope shape of all allowed EBL shapes. In most cases a single 
EBL shape represents the highest allowed shape only for a small 
wavelength interval. Consequently the envelope shape consists 
of a number of segments from different EBL shapes. One excep- 
tion occurs, if the maximum shape tested in the scan is allowed. 
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Fig. 8. H1426+428 individual results (F max > 1.5). Upper Panel: 
Two examples for the highest allowed shapes (A - solid line; B - 
dashed line) overlayed on EBL measurements and the minimum 
and maximum shapes of the scan (grey). Middle Panel: Intrinsic 
spectrum for shape A and fit functions. Both PL (solid line) and 
BPL (dashed line) are good fits; = 0.23 + 1.18 stat + 0.2 sys = 
1.61, and the BPL fit is preferred over the PL fit (P(bpl vs. pl) = 
°- 97 ); r f ; max = 2.01 + 1.94 stat +0.2 sys = 4.15andr^ L ax = -4.59+ 
6.67 stat +0.2 sys = 2.28. Lower Panel: Intrinsic spectrum for shape 
B and fit functions. Both PL (solid line) and BPL (dashed line) 
are good fits; = 1.18 + 0.24 stat + 0.2 sys = 1.63 and the BPL 
fit is not preferred over the PL fit (P(bpl vs. pl) = 0.54); 



In this case the spectrum does not constrain the EBL density. In 
general, a spectrum constrains the EBL density when the enve- 
lope shape lies below the maximum shape tested in the scan. 

Given the energy range of the spectra, the limit is only valid 
for wavelengths Au m - 
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Fig. 9. Limits on the EBL density from 1ES 1101-232 (r max > 
1.5). Grey solid curves are the minimum and maximum shapes 
of the scan; grey dashed curves are all the highest allowed EBL 
shapes for 1ES 1101-232; the thick colored curve is the corre- 
sponding envelope shape. The grey shaded area marks the wave- 
length region, in which the envelope shape is constraining the 
EBL density. 



where A mdx (E m i„) is the wavelength for which the cross section 
for pair production with the lowest energy point of the VHE 
spectrum E m i„ is maximized (following Eqn. (0). The variable 
/t thresh (£' maA ) is the threshold wavelength for pair production with 
the highest energy point of the VHE spectrum E max . This wave- 
length range roughly reflects the sensitivity of the exclusion cri- 
teria. There is, of course, some freedom of choice for the wave- 
length range, given the complexity of the criteria. 

In the following three sections the results for the individual 
prototype spectra for the realistic scan with r max > 1.5 are sum- 
marized. The results for the extreme scan with r max > 2/3 for all 
prototype spectra are presented in Sect. 15.41 

5.1. Nearby and well-measured: Mkn 501 (HEGRA) 

The envelope shape derived for the Mkn 501 spectrum is shown 
in Fig. [5] in comparison with the maximum and minimum EBL 
shapes tested in the scan. Although 7 766 674 out of 8 064 000 
EBL shapes (96.3%) can be excluded, the effective limit is only 
constraining in the 20 to 80 pm wavelength region, where it lies 
below the maximum tested shape. Note that our method is not 
only testing the overall level of the EBL density but is also sensi- 
tive to its structures. Despite the fact that, for the Mkn 501 spec- 
trum, almost all of the tested EBL shapes are excluded, certain 
types of shapes are allowed, independent of their respective EBL 
density level. This can be illustrated with an EBL shape, which 
has a power law dependency n(e) ~ e' 1 or vl v ~ AT 1 . Then the 
optical depth is independent of the energy of VHE photons and 
the intrinsic TeV blazar spec trum has the same shape as the ob- 
served one {Ah aroni anl l2001h . With such a type of shape, an al- 
lowed high EBL density level in the MIR can be constructed by 
choosing a corresponding high density level in the optical/NIR. 

The rejection power in the FIR results mainly from the hard 
intrinsic Mkn 501 spectrum above ~5TeV that often can only 
be described by an exponential or super exponential rise. Two of 
the limiting shapes, together with the resulting intrinsic spectra 
and fit functions, are shown in Fig. [6] and the corresponding fit 
results can be found in Table [3] In the table, the fit probabilities 
are given in the column "Pa". The parameters are put in paren- 
theses, if the corresponding fit has a low probability or it is not 
preferred over a fit by a function with less free parameters (result 
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Fig. 10. 1ES 1101-232 individual results (F max > 1.5). Upper 
Panel: Two examples for the highest allowed shapes (A - solid 
line; B - dashed line) overlay ed on EBL measurements and 
the minimum and maximum shapes of the scan (grey). Middle 
Panel: Intrinsic spectrum for shape A and fit functions. Both PL 
(solid line) and BPL (dashed line) are good fits; r^ x = 1.31 + 
0.13 stat + 0.1 sys = 1.53 and the BPL fit is not preferred over the 
PL fit (P(bpl vs. PL) = 0.87). Lower Panel: Intrinsic spectrum for 
shape B and fit functions. Both PL (solid line) and BPL (dashed 
line) are good fits; = 1.21 + 0.21 stat + 0.1 sys = 1.52 and the 
BPL fit is not preferred over the PL fit (P (B pl vs. pl) = 0.9488). 



of the likelihood ratio test, indicated in the column "Pref"). The 
fit with a PL function did not meet our acceptance criteria, and 
the function is omitted from the figure for the sake of legibility. 
For the two shapes displayed here, the pile up at high energies is 
already visible but not yet significant. 



10 



D. Mazin and M. Raue: New limits on the EBL density from TeV blazars 



5.2. Intermediate distance, wide energy range: H1426+428 

Using the H 1426+428 energy spectrum 5 571 772 EBL shapes 
(corresponding to 69.09% of all shapes) are excluded. The re- 
sulting envelope shape is displayed in Fig. [7] together with all 
highest allowed shapes. The limit is constraining from ~6 to 
~20 /im and the constraints are not very strong. This is illustrated 
in Fig. [8] where two representative highest allowed EBL shapes 
and the resulting intrinsic spectra are shown, as are the relevant 
fit functions. Shape A is almost as high as the maximum shape 
of the scan. Although the intrinsic spectrum is already convex, 
the relative low statistics of the spectrum even allows for an ac- 
ceptable fit by a simple PL and results in large errors on the fit 
parameters for the BPL. Consequently the EBL shape cannot 
be excluded. Shape B has a high peak in the O/IR wavelength 
region but lies below all upper limits in the FIR. The resulting 
intrinsic spectrum is best described with a simple PL with a max- 
imum slope rjna X = 1 .63 > 1.5 and is therefore allowed. 

5.3. Distant source, hard spectrum: 1ES 1101-232 

The spectrum from 1ES 1 101-232 gives the strongest constraints 
for all individual spectra, even though slightly less EBL shapes 
(7 706625, 95.57% of all shapes) are excluded than in the case 
of Mkn501. This is due to the different sensitivities in EBL 
wavelength; the 1ES 1101-232 spectrum is only sensitive up to 
10/im. The resulting maximum shapes and the envelope shape 
are shown in Fig. [9] The envelope shape constraints the EBL 
density in the wavelength range from ~0.4 to ~ 1 and clea rly 
excludes the NIR excess claimed by Matsumoto et al.l d2005l) as 
being extragalactic. At EBL wavel engths of ~2^m the lim it is 
consistent with the limit derived by Aharon ian et alj d2006al) for 
the same source with a different technique (see Sect. [6]). Two 
representative highest allowed shapes and the corresponding in- 
trinsic spectra are shown in Fig. [TOj Shape A illustrates the in- 
trinsic spectrum for a high EBL density in the UV/O, while 
shape B is the maximum shape for wavelengths around 3 fim. For 
both shapes the examined fit parameters are close to the allowed 
limits (see caption for values), as expected for highest allowed 
shapes. 

5.4. Extreme case: F max > 2/3 

The upper limits for Mkn501, H1426+428 and 1ES 1101-232 
for the extreme case with F max > 2/3 are shown in Fig. [TTlin 
comparison to the limits derived for the realistic case. In the 
case of Mkn501, a similar number of EBL shapes (94.13%) as 
in the realistic scan are excluded, but the effective limit is less 
constraining (Fig. [TTJ Upper Left Panel). For H1426+428 the 
limit almost remains at the same level, still close to the maxi- 
mum shape tested in the scan (Fig. QT| Upper Right Panel). For 
1ES 1101-232 the limit in the UV/O lies a factor of 1.2 to 1.8 
higher than the limit in the realistic case (Fig. QT| Lower Left 
Panel). In t he wavelength region fro m 2 to 4/rni, the NIR excess 
claimed by Matsum oto et al.l d2005l) is now compatible with the 
limit. In the 1 to 2fi region the limit still lies clearly below the 
claimed NIR excess. 

6. Combined results 

The number of excluded shapes for all individual spectra for the 
realistic scan are summarized in Table |4] One finds that some 
spectra (namely 1ES 2344+514, Mknl80, and PKS 2005-489) 
excluded none of the EBL shapes of the scan. This mainly is a 



Table 4. Number of excluded EBL shapes for the realistic scan 
for individual spectra (column two) as well as the number of the 
intrinsic spectra, which could not be fitted satisfactorily (column 
three). 3 



Spectrum 


#Shapes Excluded 


#Shapes No Fit 


Mkn 421 (MAGIC) 


1756 869 (21.79%) 


59 258 (0.94%) 


Mkn 421 (HEGRA) 


888 575 (11.02%) 





Mkn 421 (Whipple) 


2 287 059 (28.36%) 





Mkn 501 


7 760733 (96.24%) 


1296 (0.43%) 


1ES 2344+514 








Mkn 180 








1ES 1959+650 


1086 314(13.47%) 


423 (0.01%) 


PKS 2005-489 








PKS 2155-304 


4248 872 (52.69%) 





H 1426+428 


5 571 771 (69.09%) 


2 (0.00%) 


H 2356-309 


4 657 817 (57.76%) 





1ES 1218+304 


19 540(00.24%) 





1ES 1101-232 


7 706 624 (95.57%) 






il The percentage of these numbers to the number of all shapes 
(8 064 000) for column two and to the number of the allowed shapes for 
this spectrum for column three is also quoted. If the number of shapes 
is zero, the percentage is omitted. 




1 2 3 4 >5 



# Spectra 

Fig. 12. Number of excluded shapes as percentage of the total 
number of shapes vs. the number of spectra, that excluded the 
shape. A large fraction of shapes (98.08%) is excluded by more 
than one spectra. 



result of the low statistics in the spectra, hence the large errors 
on the fit parameters. The highest number of excluded shapes 
and the strongest constraint in the scan result from the previ- 
ously presented prototype spectra Mkn501, H 1426+428, and 
1ES 1101-232. 

We now combine the results from all spectra by excluding 
all shapes of the scan, which are excluded by at least one of the 
spectra. First we present results for the realistic scan. A large 
fraction of the shapes (98.08%) are excluded by more than one 
spectra (see Fig. [12), which strengthens our results. In total, 
we exclude 8 056718 shapes, which leaves us with only 7 282 
(0.01 %) allowed shapes. The resulting maximum shapes (dashed 
lines) and the envelope shape for the combined results are shown 
in the upper panel of Fig. Q~3]in comparison to the results from 
the individual prototype spectra. In the optical-to-near-infrared 
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Fig. 11. Result for individual spectra, comparing the realis- 
tic case with r max > 1 .5 (thick colored solid line) with the 
extreme case with r max > 2/3 (thick colored dashed line), 
overlaid on the EBL measurements and the minimum and 
maximum shapes of the scan (grey markers and lines). The 
grey boxes indicate the validity range of the limit for the 
realistic scan as in Figs. HIT] and [9] Note that the validity 
range of the limits from the extreme scan are quite similar. 
Upper Left Panel: Results for Mkn 501 (HEGRA). Upper 
Right Panel: Results for H1426+428. Lower Left Panel: 
Results for 1ES 1101-232. 



the combined limit follows exactly the limit derived from the 
1ES 1101-232 spectrum. For longer wavelengths, however, the 
combined limit lies significantly below the limits derived from 
individual spectra. This is particularly striking in the MIR, where 
the H 1426+428 spectrum provided only weak limits and the 
combined limit now gives considerable constraints. This can be 
understood from the fact that, for a high EBL density in the MIR, 
a high density in the optical to NIR is needed, so that the spectra 
(Mkn 501, H 1426+428) will not get too hard. These high den- 
sities in the optical-to-NIR are now excluded by the 1ES 1110- 
232 spectrum, which therefore result in stronger constraints in 
the MIR. 

The combined limit for the realistic scan in comparison 
to the direct measurements and limits is shown in the lower 
panel of Fig. [13] In the optical-to-NIR one finds that the 
combined limit lies signif i cantly below the claimed NIR ex- 
cess by iMatsumoto et alj d2005l). It is compa t ible with the 
detecti o ns reported by iDwek & Arendtl (1 19981). iGoriian et al.l 
d2000l) . IWright & Reesd (120001) . and ICambresv et all (l200lT 
In the same figure the limit reported by the H.E.S.S. collab- 
oration for 1ES 1101-2 32 derived with a different technique 
(Ahar onianet al . 2006a) is shown and at wavelengths around 2 - 
3 pm it is in good agreement with the limit derived here; but for 
shorter wavelengths our limit lies significantly higher. While for 
the H.E.S.S. limit comparable exclusion criteria for the spectra 
were used, a fixed reference shape scaled in the level of the EBL 
density was used to calculate the EBL limi0, which presumably 
causes the differences at smaller wavelengths. In the MIR to FIR 
our combined limit lies below all previously reported upper lim- 
its from direct measurements and fluctuation analysis. For EBL 
wavelengths greater than 2 pm, it is only about a factor of 2 to 2.5 
higher than the absolute lower limit from source counts, leaving 
very little room for additional contributions to the EBL in this 



3 In addition to the fixed shape, several oth er possible EBL compo - 
nents, like a bump in the UV, were examined ( Aharonian et al. 2006a). 



wavelength region. In the FIR it lies more than a facto r of ~ 2 
below the claimed detection by Finkb einer et al.l d2000t) . 

In the realistic scan, after combining the results from 
Mkn 501, H 1426+428, and 1ES 1110-232, adding more spec- 
tra does not further strengthen the limit (though marginal more 
shapes are excluded). Hence, for the extreme scan, we only com- 
bine the results from these three spectra. The limit for the ex- 
treme scan is shown in Fig. [14] in comparison to the limit from 
the realisistic scan. The limit from the extreme scan lies for 
all EBL wavelengths above the limits from the realistic scan 
and is a factor of 2.5 to 3.5 higher than the lower limits from 
source counts. In the optical-to-NIR the combined limit again 
follows the limit derived from the 1ES 1101-232 spectrum. As 
stated for the ind i vidua l spectrum the NIR excess claimed by 
IMatsumoto et alj d2005l) for wavelengths above 2 pm is now 
compatible with the limit, but for shorter wavelengths it is still 
clearly excluded. 

In Fig. [15] the limits derived in this paper are shown in 
compa r ison to predic t ions f rom EBL models by IKneiske et al.l 
d2002l) . IStecker et alj d2006l) . and iPrimack et alf d2005l) (at z = 
0). Most models lie below o ur realisitc limits. Only the high 
model of IKneiske et al.l d2002 ) is slightly h i gher i n the NIR, and 
the fast evolution model of IStecker et alj (120061) is marginally 
higher in the MIR. These two models are also excluded when 
applying the criteria from Sect. [4] while the other m odels are al- 
lowed. Not eworthy is that the Primack et al.l d2005l) and the low 
model from IKneiske et all d2002l) are below the lower limits de- 
rived from galaxy counts in the MIR and FIR. 



7. Summary and conclusions 

In this paper we present a new method of constraining the den- 
sity of the extragalactic background light (EBL) in the optical- 
to-far-infrared wavelength regime using VHE spectra of TeV 
blazars. We derive strong upper limits, which are only a factor 2 
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Fig. 13. Combined results for 
the realistic scan. Upper Panel: 
All highest allowed shapes of 
the combined scan (dashed grey 
lines) and the corresponding en- 
velope shape (solid black line) in 
comparison to the limits for in- 
dividual spectra: Mkn501 (solid 
blue line), H 1426+428 (solid 
yellow line), and 1ES 1101-232 
(solid red line). The minimum 
and maximum shapes of the 
scan are also shown (grey lines). 
Lower Panel: The combined 
limit from the realistic scan 
(solid black line) in compari- 
son to direct measurements and 
limits (grey marker). The grey 
dashed curve ar ound 2 //m is the 
limit de rived by 1 Aharonian et af] 
(20063). 
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Fig. 14. Combined results from 
the extreme scan (dashed black 
line) in comparison to the re- 
sult from the realistic scan (solid 
black line). Grey lines are the 
minimum and the maximum 
shapes tested in the scan. Grey 
markers are direct measurements 
and limits. 



to 3 higher than the lower limits determined from source counts. 
Unlike what has been done in many previous studies of this kind, 
we do not use one or a few pre-defined EBL shapes or models, 
but rather a scan on a grid in EBL density vs. wavelength. This 
grid covers wavelengths from 0.1 pm to 1000//m and spans in 
EBL density from the lower limits from source counts to the up- 
per limits determined from direct measurements and fluctuation 



analyses. By iterating over all scan points, we test a large set of 
different EBL shapes (8 064 000 in total). Each EBL shape is de- 
scribed by a spline function, which allows us to calculate the op- 
tical depth of VHE y-rays for this shape via a simple summation 
instead of solving three integrals numerically. The resolution of 
the scan for sharp bumps or dips in the EBL is limited by the 
choice of the grid and the order of the spline. 
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Fig. 15. EBL limits from the re- 
alistic (thick solid black curve) 
and the extreme scan (thick 
dashed black curve) in compar- 
ison to different EBL models at 
z = (blue cur ves): updated 
versio n of the iKneiske et al.l 
(120021) high and low mo dels 
(solid), IStecker et al.l (120061) fast 
and baseli ne evolution models 
(dashed), iPrimack et ail d2005l) 
(dashed dotted). 
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To test these EBL shapes, we use the measured VHE y-ray 
spectra from distant sources, which carry an imprint of the EBL 
from attenuation via pair production. We included spectra from 
all known TeV blazars (13 in total), making this the most com- 
plete study of this type up to now. All spectra are scrutinized 
with the same robust algorithm. For each spectrum and each 
EBL shape, the intrinsic spectrum is calculated and an analyt- 
ical description of the intrinsic spectrum is determined by fitting 
several functions. The fit parameters and errors of the functions 
are subsequently used to evaluate whether the intrinsic spectrum 
is physically feasible. We use two conservative criteria from the 
theory of VHE y-ray emission in TeV blazars. The EBL shape 
is excluded if (i) a part of the intrinsic spectrum is harder than 
a theoretical limit or (ii) the intrinsic spectrum shows a signifi- 
cant pile-up at high energies. Since there is some spread in the 
predictions from theory, we included two independent scans for 
two theoretical limits on the hardness of the spectrum: (1) the 
realistic case for a photon index of T > 1.5 and (2) the extreme 
case with F > 2/3. 

We present limits derived from individual VHE y-ray spec- 
tra. The strongest constraints result from the 1ES 1 101-232 spec- 
trum in the wavelength range from 0.8 and 10 /mi, mainly due to 
the hard spectrum and its large distance to earth. An extragalac- 
tic origin of the claimed excess in the NIR between 1 and 2 /mi 
( Matsumo to et al.l 120051) can be excluded even by the extreme 
sc an. This result confirms th e results from the recent publication 
of lAharonian et all (l2006al) . In the FIR (20 and 80/mi), strong 
upper limits are provided by the nearby TeV blazar Mkn501. 
H 1426+428, situated at an intermediate distance, provides some 
constraints on the EBL density in the MIR (between 1 and 
15 /mi), connecting the upper limits derived from 1ES 1 101-232 
and Mkn501. Most of the other spectra provide constraints in 
certain wavelength ranges, but they are not as strong as the lim- 
its from the previously described sources. 

By combining the results from all spectra, we find that the 
upper limits become much stronger compared to the limits de- 
rived for individual spectra, especially in the wavelength range 
between 4 and 60 /mi. In the wavelength region from 2 to 80 /mi, 
the combined limit for the realistic case lies below the upper lim- 
its derived from fluctuation an alyses of the direct measurements 
dKashlinskv & Odenwaldll2000l) and is just a factor of 2 to 2.5 
above the absolute lower limits from source counts. This makes 
it the most constraining limit in the MIR to FIR region so far. As 
expected, the upper limits from the extreme scan are less con- 



straining (factor of 2.5 to 3.5 above the lower limit from source 
counts), but still substantial over this wide wavelength range. 
The derived upper limits can be interpreted in two different 

ways: 

1 . The EBL density in the optical-to-FIR is significantly lower 
than suggested by direct measurements, and the actual EBL 
level seems to be close to the existing lower limits. It can thus 
be concluded that experiments like HST, ISO, and Spitzer re- 
solve most of the sources in the universe. This would indicate 
that there is very little room left for a significant contribution 
of heavy and bright Population III stars to the EBL density 
in this wavelength region. 

2. The assumptions used for this study are not correct. This 
would require a revision of the current understanding of TeV 
blazar physics and models, which has so far been fairly suc- 
cessful in modeling multi-wavelength data from the radio to 
the VHE for all detected sources. 

We estimate that the main contribution to the systematic er- 
rors arises from the minimum width of the EBL structures that 
can be resolved by the scan. Although the choice of the grid 
spacing, which defines the minimum width of the EBL struc- 
tures, is physically well motivated, there are arguments that even 
thinner structures could be realized in nature due to, say absorp- 
tion effects. We tested a 20% smaller grid spacing, resulting in 
much less realistic, but still possible, EBL structures, with the 
1ES 1 101-232 spectrum. The resulting limits on the EBL are 20 
to 30% higher than the ones presented in Section 5.3. We there- 
fore conclude that the overall contribution from the grid setup to 
the systematic error is at most 30%. 

Another contribution to the systematic error originates from 
not considering the evolution of the EBL in our method. To es- 
timate the effect of the evolution of the EBL for the most dis- 
tant sources, we calculate the late contribution to the EBL in 
the redshift interval bet ween z = and z = 0.2. Here, we uti- 
lize the EBL model by Knei ske et al.l d2002l) (updated version 
is used, Kneiske private communication). We obtain a wave- 
length dependent contribution, which has a maximum value of 3- 
4nWrrT 2 sr~ 1 in the optical to NIR (the relevant wavelength re- 
gion for the distant sources in our scan). Given that these values 
are derived with the extreme assumption that all late emission 
occured instantaneously at z=0, we estimate the systematic er- 
ror arising from the EBL evolution to 10% of the derived limits. 
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Note that Aharonia n et al.l d2006a) estimated an error of <10% 
for the most distant TeV blazar in our study 1ES 1 101-232. 

To study the effect of the uncertainties in the absolute en- 
ergy scale of IACT measurements, we calculate a limit on the 
EBL utilizing the energy spectra from Mkn501, H 1426+428, 
and 1ES 1 101-232 shifted by 15% to lower energies (for F > 1.5, 
realistic scan). Since the optical depth is increasing with energy, 
the energy shift reduces the attenuation of the spectra, especially 
reducing the effect of a possible pile-up at the highest energies. 
The difference between the limit derived from the energy shifted 
spectra and the realistic limit is negligible in the optical to NIR 
(<6/mi), while there is a moderate effect in the MIR to FIR on 
the level of 1 -4.5 nWm^sf 1 . We therefore conclude that the 
systematic error in the MIR to FIR (A > 6 ^m) introduced by the 
uncertainties in the absolute energy scale of IACT measurements 
is close to 10-45% of our realistic limit. 

Adding the individual errors quadratic ally, we finally con- 
clude that the (very conservative) systematic error on the upper 
limit is about 32% in the optical-to-NIR and 33-55% in MIR to 
FIR. 

The strong constraints derived in this paper only allow for 
a low level of the EBL in the optical to the far-infrared, sug- 
gesting that the universe is more transparent to VHE y-rays than 
previously thought. Hence, we expect detections of many new 
extragalactic VHE sources in the next few years. Further multi- 
wavelength studies of TeV blazars will improve the understand- 
ing of the underlying physics, which will help to refine the ex- 
clusion criteria for the VHE spectra in this kind of study. The 
upcoming GLAST satellite experiment, operating in an energy 
range from 0.1 up to ~100 GeV, will allow such studies of the 
EBL to be extended to the ultraviolet-to-optical wavelength re- 
gion. A detailed analysis of the impact of our constraints on the 
contribution from Population III stars to the EBL is in prepara- 
tion. 

Acknowledgements. The authors would like to thank J. Ripken for fruitful 
discussion of splines and the careful reading of the manuscript. The authors 
also wish to thank M. Beilicke, R. Cornils, F. Goebel, N. Goetting, and 
G. Heinzelmann for stimulating comments and help in preparing this manuscript. 
The authors also would like to thank the anonymous referee for the helpful com- 
ments, which improved the manuscript. D. M. acknowledges the financial sup- 
port of the MPI fur Physik, Munich. M. R. acknowledges the financial support 
of the University of Hamburg and the BMBF. This research made use of NASA's 
Astrophysics Data System. 

Appendix A: Likelihood Ratio Test 

The likelihood ratio test (lEadie et alJll988h is a standard statis- 
tical tool to test between two hypotheses whether an improve- 
ment to a fit quality (quantified by corresponding x 2 values) is 
expected from a normal distribution or if it is significant. By fit- 
ting two functional forms to the intrinsic spectrum, one obtains 
values of the likelihood functions La and Lb- If hypothesis A 
is true, the likelihood ratio R = -\\\(LaI Lb) is approximately 
X 2 distributed with N degrees of freedom. N is the difference 
between numbers of degrees of freedom of hypothesis A and 
hypothesis B. One defines a probability 

P= J p(*W (A.i) 



where p (x 2 ) is the x 2 probability density function and R meas the 
measured value of R. Hypothesis A will be rejected (and hypoth- 
esis B will be accepted) if P is greater than the confidence level, 
which is set to 95%. 
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Table 2. Analytical functions that are used to fit the intrinsic TeV blazar spectrum. 



# Description 



Abbreviation Formula /(£) = dN/dE 



Parameters to evaluate 



2 pfL 



1 simple power law 

2 broken power law with transition region 

3 broken power law with transition region 
and super-exponential pile-up 

4 double broken power law with transi- 
tion regions 

5 double broken power law with transi- 
tion regions and super-exponential pile- 
up 



PL 
BPL 



DBPL 



BPLSE NoE^[l + (0] f exp(^) 



DBPLSE DBPL x exp ) 



BPL 

2 



X 



pDBPL pDBPL pDBPL 



x 1 



